import struct
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_mldata
from sklearn.base import BaseEstimator
from sklearn.linear_model import SGDClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score,cross_val_predict
from sklearn.multiclass import OneVsOneClassifier
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import KNeighborsClassifier
from sklearn.metrics import confusion_matrix,precision_score,recall_score,f1_score,precision_recall_curve,\
    roc_curve,roc_auc_score

with open('./MNIST_data/train-images-idx3-ubyte','rb') as f:
    buffer = f.read(4*4) #4个int
    # 解包
    head = struct.unpack('>iiii',buffer)
    # print(head)
    length = head[1]*head[2]*head[3]
    # print(length)
    buffer = f.read(length)
    data = struct.unpack('>{}B'.format(length),buffer)
    # print(data)
    # images = np.reshape(data,(head[1],head[2],head[3]))
    # print(images)
    # for i in range(5):
        # plt.imshow(images[i],cmap='Accent')
        # plt.show()
    mnist = fetch_mldata('MNIST original',data_home='./')
    # print(np.array(mnist['data']).shape)
    x,y=mnist['data'],mnist['target']
    some_digit = x[36000]
    some_digit_image = some_digit.reshape(28,28)
    # plt.imshow(some_digit_image,cmap=matplotlib.cm.binary)
    # plt.show()
    x_train,x_test,y_train,y_test = x[:60000],x[60000:],y[:60000],y[60000:]
    shuffle_index = np.random.permutation(60000)
    x_train,y_train = x_train[shuffle_index],y_train[shuffle_index]

    # 训练一个二元分类器
    # 识别数字5，二元分类5或者非5
    # 创建目标向量
    y_train_5 = (y_train == 5)
    y_test_5 = (y_test == 5)

    # # SGD梯度下降分类器，适合非常大的数据集，独立处理训练集数据，一次一个，适合在线学习
    sgd_clf = SGDClassifier(random_state=42)
    sgd_clf.fit(x_train,y_train_5)
    # print(sgd_clf.predict([some_digit]))
    # 评估分类器比评估回归器要困难得多
    # score = cross_val_score(sgd_clf,x_train,y_train_5,cv=3,scoring='accuracy')
    # print(score)
    # class NeverClassifier(BaseEstimator):
    #     def fit(self,X,y=None):
    #         pass
    #     def predict(self,X):
    #         return np.zeros((len(X),1),dtype=bool)

    # never_5_clf = NeverClassifier()
    # print(cross_val_score(never_5_clf,x_train,y_train_5,cv=3,scoring='accuracy'))
    # 混淆矩阵
    # y_train_pred = cross_val_predict(sgd_clf,x_train,y_train_5,cv=3)
    # print(confusion_matrix(y_train_5,y_train_pred))
    #PR曲线
    #第一行一列真负列，第一行二列假正列，第二行第一列假负列，第二行第二列真正列
    #正类预测的准确率被称为分类器的精度 真正/真正+假正
    #真类预测的准确率 召回率 真正/真正+假负
    # print(precision_score(y_train_5,y_train_pred))
    # print(recall_score(y_train_5,y_train_pred))
    #F1谐波平均值
    #F1分数对那些具有相近精度和召回率，分类器更有利，这不一定符合你的期望
    # print(f1_score(y_train_5,y_train_pred))
    #训练一个分类器检测儿童可以放心观看的视频，你可能要求拦截了很多的视频，低召回率，保留下来的都是安全的视频，高精度
    #不能同时增加精度并减少召回率，反之亦然
    #阀值左侧为负类，右侧为正类，提高阀值精度增高，召回率降低，skl不能直接设置阀值，可以访问决策分数，SGDClassifier默认阀值为0
    #返回决策值
    # y_scores = sgd_clf.decision_function([some_digit])
    # print(y_scores)
    #设置阀值
    threshold = 0
    # y_some_digit_pred = (y_scores>threshold)
    #返回决策分数，而不是预测结果
    # y_scores = cross_val_predict(sgd_clf,x_train,y_train_5,cv=3,method='decision_function')
    # precisions,recalls,thresholds = precision_recall_curve(y_train_5,y_scores)
    # def plot_precision_recall_vs_threshold(precisions,recalls,thresholds):
    #     plt.plot(thresholds,precisions[:-1],'b--',label = 'Precision',linewidth=2)
    #     plt.plot(thresholds, recalls[:-1], 'g-', label='Recall', linewidth=2)
    #     plt.xlabel('Treshold',fontsize=16)
    #     plt.legend(loc='upper left',fontsize=16)
    #
    # plt.figure(figsize=(8,4))
    # plot_precision_recall_vs_threshold(precisions,recalls,thresholds)
    # plt.xlim([-700000,700000])
    # plt.show()
    # y_train_pred_90 = (y_scores>70000)
    # print(precision_score(y_train_5,y_train_pred_90))
    # print(recall_score(y_train_5,y_train_pred_90))
    # plt.rcParams['font.sans-serif']=['SimHei']
    # def plot_precision_vs_recall(precisions,recalls):
    #     plt.plot(recalls,precisions,'b-',linewidth=2)
    #     plt.xlabel('召回',fontsize=16)
    #     plt.ylabel('精度',fontsize=16)
    #     plt.axis([0,1,0,1])
    # plt.figure(figsize=(8,6))
    # plot_precision_vs_recall(precisions,recalls)
    # plt.show()
    #总结 ：获得了一个90%精度的分类器，但如果召回太低，精度再高，也不怎么有用
    #ROC曲线
    #本质是真正类率tpc和假正类fpr（错误的分为正类的负类实力比例）
    #与召回/精度曲线非常相似
    # fpr,tpc,threshold = roc_curve(y_train_5,y_scores)
    # plt.rcParams['font.sans-serif']=['SimHei']
    #
    #
    # def plot_roc_curve(fpr,tpc,label=None):
    #     plt.plot(fpr,tpc,linewidth=2,label=label)
    #     plt.plot([0,1],[0,1],'k--')
    #     plt.axis([0,1,0,1])
    #     plt.xlabel('假正类率',fontsize=16)
    #     plt.ylabel('真正类率',fontsize=16)
    # plt.figure(figsize=(8,6))
    # plot_roc_curve(fpr,tpc)
    # plt.show()
    #计算曲线下面积AUC
    # print(roc_auc_score(y_train_5,y_scores))
    #召回率越高，分类器的假正类FPR就越多
    #虚线表示纯随机分类的ROC曲线，好的分类器应该远离这条线，向左上角
    #是选择PR曲线还是ROC曲线，正类如果非常少或更关注假正而不是假负类选择PR，反之ROC
    #训练随机森林分类器，比较SGD分类器的ROC曲线和ROC AUC分数
    #获取每个实例的分数
    #RandomForestClassifier没有decisi_functi（），但是拥有dict_proda（）方法，skl中分类器都有这两个中的一个
    #dict_proda返回一个矩阵，每行一个实例，每列代表一个类别的概率，比如这个图片70%是5
    # forest_clf = RandomForestClassifier(n_estimators=10,random_state=42)
    # y_prooba_forest = cross_val_predict(forest_clf,x_train,y_train_5,cv=3,method='predict_proba')
    #绘制ROC曲线需要决策值不是概率，直接使用正类的概率作为决策值
    # y_scores_forest = y_prooba_forest[:,1]
    # fpr_forest,tpc_forest,thresholds_forest = roc_curve(y_train_5,y_scores_forest)
    # plt.figure(figsize=(8,6))
    # plt.plot(fpr,tpc,'b:',linewidth=2,label='SGD')
    # plt.plot(fpr_forest,tpc_forest,linewidth=2,label='Random Forest')
    # plt.xlabel('假正类率', fontsize=16)
    # plt.ylabel('真正类率',fontsize=16)
    # plt.legend(loc='lower right',fontsize=16)
    # plt.show()
    #Rand 比SGD好很多，ROC AUC的分数也很高很多
    # print(roc_auc_score(y_train_5,y_scores_forest))
    # y_train_pred_forest = cross_val_predict(forest_clf,x_train,y_train_5,cv=3)
    #精度和召回率也很高
    # print(precision_score(y_train_5,y_train_pred_forest))
    # print(recall_score(y_train_5,y_train_pred_forest))
    #总结 选择合适的指标利用交叉验证来对分类进行评估
    #选择满足需求的精度/召回率权衡
    #使用ROC曲线和ROC AUC分数比较多个模型
    #多类别分类器
    #尝试5之外的检测
    #多类别分类器区分两个以上的类别
    #随机森林和朴素贝叶斯可以直接处理多个类别
    #支持向量机和线性分类器只可以处理二元分类器
    #解决方案
    #将图片分类0到9，训练10个二元分类器，每个数字一个，检测一张图片时，获取每个分类器的决策分是，那个最高属于哪个，称为一对多ova
    #为每一对数字训练一个二元分类器，区分0，1区分0，2区分1，2称为一对一ovo策略，存在n个类型，需要N*（n-1）/2个分类器，最后看哪个类别获胜最多
    #优缺点
    #ovo只需要用到部分训练集对其必须区分两个类型进行训练
    #对于较小训练集ovo比较有优势，大训练集ova更常用，ova速度快
    #skl检查到使用二元分类算法进行多类别分类任务，会自动运行ova，svm分类器除外
    # sgd_clf.fit(x_train,y_train)
    # sgd_clf.predict([some_digit])
    #内部实际上训练了10个二元分类器，获得图片的决策分数，然后选择了分数最高的类别
    #返回10个分数，每个类别一个
    # some_digit_scores = sgd_clf.decision_function([some_digit])
    # print(some_digit_scores)
    # sgd_clf.classes_[np.argmax(some_digit_scores)]
    #使用ovo策略，一对一或者一对多
    # ovo_clf = OneVsOneClassifier(SGDClassifier(max_iter=5,tol=np.infty,random_state=42))
    # ovo_clf.fit(x_train,y_train)
    # ovo_clf.predict([some_digit])
    # print(len(ovo_clf.estimators_))
    #使用随机森林
    # forest_clf.fit(x_train,y_train)
    # forest_clf.predict([some_digit])
    #随机森林直接将实例分为多个实例，调用predict_proba（）可以获得分类器将每个实例分为每个类别的概率列表
    # print(forest_clf.predict_proba([some_digit]))
    #评估分类器
    #使用交叉验证评估SGD的准确率
    # print(cross_val_score(sgd_clf,x_train,y_train,cv=3,scoring='accuracy'))
    #将输入进行简单缩放，可以得到准确率90%以上
    scaler = StandardScaler()
    x_train_scaled = scaler.fit_transform(x_train.astype(np.float64))
    # print(cross_val_score(sgd_clf,x_train_scaled,y_train,cv=3,scoring='accuracy'))
    #错误分析
    #项目流程
    #1.探索数据准备的选项
    #2.尝试多个模型
    #3.选择最佳模型并用GridSearchCV对参数进行预测
    #4.尽可能自动化
    #确定了一个相对合适的模型，进一步优化，分析其错误类型
    #1.查看混淆矩阵
    #2.使用cross_val_predict()进行预测
    #3.使用confusion_matrix()
    y_train_pred = cross_val_predict(sgd_clf,x_train_scaled,y_train,cv=3)
    conf_mx = confusion_matrix(y_train,y_train_pred)
    # print(conf_mx)
    #使用matplotlib的matshow函数来查看混淆矩阵的图像表示
    # plt.matshow(conf_mx,cmap=plt.cm.gray)
    # plt.show()
    #看起来不错，大多数图片都在主对角线上，说明它们被正确分类
    #数字5看起来比较暗，说明1，数字5图片比较少 2.分类器在数字5上执行效果不如其他数字好
    #假设把焦点放在错误上，为取得错误率，而不是错误绝对值，需要将混淆矩阵中每个值除以相应类别中的图片数量
    # row_sums = conf_mx.sum(axis = 1,keepdims=True)
    # norm_conf_mx = conf_mx/row_sums
    #用0填充对角线，只保留错误，重新绘制
    # np.fill_diagonal(norm_conf_mx,0)
    # plt.matshow(norm_conf_mx,cmap=plt.cm.gray)
    # plt.show()
    #每行代表实际类别，没列代表预测类别
    #8，9列比较亮，说明许多图片被错误的分类为数字8，9
    #8，9行也偏亮，说明数字8和9经常跟其他数字混淆
    #有些暗，比如行1，大多数数字1都被正确的分类，一些和8混淆
    #3和5是错误最多的
    #结论
    #改进数字8和9的分类
    #修正数字3和5的混淆
    #如何优化分类器
    #尝试多收集这些数字的训练集
    #开发一些新特征来改进分类器
    #优化分类器算法
    #使用pillow和opencv对图片预处理，让显示模型更突出
    #分析单个错误
    # def plot_digits(instances,image_per_row=10,**options):
    #     size=28
    #     image_per_row=min(len(instances),image_per_row)
    #     images=[instance.reshape(size,size) for instance in instances]
    #     n_rows = (len(instances)-1)//image_per_row+1
    #     row_images=[]
    #     n_empty = n_rows*image_per_row-len(instances)
    #     images.append(np.zeros((size,size*n_empty)))
    #     for row in range(n_rows):
    #         rimages = images[row*image_per_row : (row+1)*image_per_row]
    #         row_images.append(np.concatenate(rimages,axis=1))
    #     image = np.concatenate(row_images,axis=0)
    #     plt.imshow(image,cmap=matplotlib.cm.binary,**options)
    #     plt.axis('off')
    # #查看数字3和5的例子
    # cl_a,cl_b=3,5
    # x_aa=x_train[(y_train==cl_a)&(y_train_pred==cl_a)]
    # x_ab=x_train[(y_train==cl_a)&(y_train_pred==cl_b)]
    # x_ba=x_train[(y_train==cl_b)&(y_train_pred==cl_a)]
    # x_bb=x_train[(y_train==cl_b)&(y_train_pred==cl_b)]

    # plt.figure(figsize=(8,8))
    # plt.subplot(221);
    # plot_digits(x_aa[:25],image_per_row=5)
    # plt.subplot(222);
    # plot_digits(x_ab[:25], image_per_row=5)
    # plt.subplot(223);
    # plot_digits(x_ba[:25], image_per_row=5)
    # plt.subplot(224);
    # plot_digits(x_bb[:25], image_per_row=5)
    # plt.show()
    #左侧两个被分类为3的图像
    #右侧两个被分类为5的图像
    #大多数错误被分类的图片看起来还是非常明显的错误
    #原因：SGD是一个线性模型，它所做就是为每个像素分配一个各个类别的权重，当它看到新的图像，将加权后的像素强度汇总，从而得到一个分数进行分类
    #数字3和5在一部分像素位上有区别，所以分类器很容易将其弄混
    #通过上面图像，如果书写3的连接点左移，分类器可能将其分类为数字5，这个分类器对图像位移和旋转敏感
    #减少混淆的方法之一，就是对图像进行预处理，确保位于中心位置并且没有旋转
    #多标签分类
    #为每个实例产生多个类别，例如：照片识别多个人脸
    #分类器经过训练可以识别小红，小白，小军，一张照片里有小红，小白
    #经过分类器，应该输出【1，1，0】是小红，是小白，不是小军
    #输出多个二元标签的分类系统称为多标签分类系统
    # y_train_large = (y_train>=7)
    # y_train_odd = (y_train % 2 ==1)
    # y_multilabel = np.c_[y_train_large,y_train_odd]
    knn_clf = KNeighborsClassifier()
    # knn_clf.fit(x_train,y_multilabel)
    #knn支持多标签分类，不是所有分类器都支持多标签分类
    # print(knn_clf.predict([some_digit]))
    #评估多标签分类器方法很多，方法之一就是测量每个标签的F1分数，或者其他二元分类器指标，然后简单平均
    # y_train_knn_pred=cross_val_predict(knn_clf,x_train,y_multilabel)
    # print(f1_score(y_multilabel,y_train_knn_pred,average='macro'))
    #上面假设了所有标签都等同重要，也可以给每个标签设置一个权重（该目标标签实例的数量），设置avera=‘weighted’
    #多输出分类
    #例子构建一个系统去除图片中的噪声，输入一张有噪声的图片，他将输入一张干净的数字图片，分类器输出是多个标签，一个像素一个标签，每个标签多个值0到255
    #增加噪声，目标将图片还原为原始图片，创建训练集和测试集
    # noise = np.random.randint(0,100,(len(x_train),784))
    # x_train_mod = x_train + noise
    # noise = np.random.randint(0,100,(len(x_test),784))
    # x_test_mod = x_test+noise
    # y_train_mod = x_train
    # y_test_mod = x_test
    #
    # some_index = 5500
    # # plt.subplot(121);
    # # plt.imshow(x_test_mod[some_index].reshape(28,28),cmap=matplotlib.cm.binary)
    # # plt.subplot(122);
    # # plt.imshow(y_test_mod[some_index].reshape(28, 28), cmap=matplotlib.cm.binary)
    # # plt.show()
    # knn_clf.fit(x_train_mod,y_train_mod)
    # clean_digit = knn_clf.predict([x_test_mod[some_index]])
    # plt.imshow(clean_digit.reshape(28,28),cmap=matplotlib.cm.binary)
    # plt.show()

















